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Abstract 

The solvation of Al^+ and its hydrolyzed species in water clusters has been studied 
by means of ab initio molecular dynamics simulations. The hexa-hydrate Al(H20)g^ 
ion formed a stable complex in the finite temperature cluster simulation of one alu- 
minum ion and 16 waters. The average dipole moment of strongly polarized hydrated 
water molecules in the first solvation shell of A1(H20)q^ was found to be 5.02 De- 
bye. The deprotonated Al(H20)2(OH)j complex evolves into a tetra-coordinated 
A1(0H)J aluminate ion with two water molecules in the second solvation shell 
forming hydrogen bonds to the hydroxyl groups in agreement with the observed 
coordination. At high temperature and for low water coordination protons in the 
first solvation shell are very mobile leading to the formation of hydrolysis species 
consistent with the acidity of Al^+ in water. 



1 Motivation 



The aqueous chemistry of aluminum, the most abundant metal in the earth's 
crust, remains a subject of fundamental research after many years of inten- 
sive research[l]. Due to the increased acidity of many natural waters and the 
corresponding increase in solubility of aluminum containing minerals, environ- 
mental and health issues have further focus attention on the solution behavior 
of this toxic material [2,3]. In addition aluminum hydrolysis products have 
many practical applications, ranging from the pharmaceutical design to pu- 
rification of water [1]. 
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The aluminum ion, AF"*", exists as an octahedral (6- fold coordinated) hexahy- 
drate ion, A1(H20)5^ in acidic solutions [4]. Over a narrow pH range, 5.5 < 
pH < 6.2, the hexahydrate undergoes hydrolysis by four successive deprotona- 
tions, producing species with uniformly decreasing coordination numbers[5]. 
The final product of the hydrolysis process, the aluminate ion, A1(0H)J, is 
tetrahedrally (4-fold) coordinated [5]. The coordination behavior of ions in so- 
lution is a sensitive function of the properties of the ions and the solvating 
molecules and varies widely from species to species. For example, the hy- 
drolyzed aqueous ferric ion, Fe^"*", remains six-coordinated after four deproto- 
nations [5]. The change in coordination of A1(H20)2(0H)J from octahedral 
to tetrahedral form occurs to accommodate changes in ligands charges and 
polarizations [5] . Because of the small relative size of the AF"*" ion, this change 
in coordination both strongly polarizes the OH bond of water ligands and 
stabilizes the formation of 0H~. Both these effects should increase acidity 
constant of Al^+ in water solutions, which is about 10"^ at room tempera- 
ture [6]. However, the microscopic description of these processes, especially 
the mechanisms of hydrolysis, remains incomplete. 

We will show in this Letter using ab initio molecular dynamics (AIMD) [7] that 
this strong coordination interaction results in highly inhomogeneous electron 

charge distribution in the first solvation shell of solvating water molecules. D. 
Marx et al. have successfully applied a similar method to the studies of the 
structural and dynamical properties of bivalent aqueous Be^+ ion [8]. However, 
there have been a number of efforts to simulate highly charged ions in solu- 
tion using traditional molecular dynamics methods. Given the complexity of 
many-body interactions in such systems with highly charged ions and hydro- 
gen bonded water molecules, the development of phenomenological potentials 
describing the interaction of trivalent ions (Al^+, Cr^+) with water and further 
classical molecular dynamics simulations often neglect hydrolysis [9,10]. If the 
goal is to model the solvation of Al^+ in the practically important range of 
pH = 4 — 7, hydrolysis has to be taken into account to correctly describe the 
stability of the surrounding aluminum water network [11,12]. 



2 Technical details 

In this Letter we illustrate the essential features of AF+ hydrolysis in terms 
of the hydration numbers and structural parameters of the species in aqueous 
solutions, namely, Al(H20)g+ and Al(OH)j. We will show that the efficient 
local density approximation (LDA), which has been used throughout the pa- 
per, provides accurate structural information for both species (see Section 3). 
Vosko et aFs parameterized form for the exchange-correlation energy has been 
used to implement LDA [13]. Since the potential energy surfaces of small clus- 
ters usually have many local minima, we have used the method of dynamical 
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simulated annealing to find the optimal structures on the ground state poten- 
tial surface. Such a search for a global minima would be difficult to do using 
only conventional total energy methods such as steepest descent minimization. 
To get accurate energetics, the Perdew-Burke-Ernzerhof 1996 (PBE96) gen- 
eralized gradient approximation [14] has been employed in a postprocessing 
mode. 

The valence electronic wave functions were expanded in plane waves, and their 
interactions with nuclei and the core electrons were described through gener- 
alized norm-conserving Hammann pseudopotentials [15]. The nonlocal part of 
the pseudopotentials was modified to a completely separable form as suggested 
by Kleinman and Bylander [16]. Since the original Hammann pseudopotential 
requires too high cut-off energy for oxygen, a softer potential was constructed 
by increasing the core radii. Using this potential, the equilibrium bond dis- 
tances and the binding energies of single water molecule and the water dimer 
were tested against known LDA results [17] at the cutoff energy Ecut = 100 
Ryd (see Table 1). Good agreements in these tests indicated that our soft 
pseudopotential for the oxygen is accurate. This choice of the cut-off energy 
required about 130,000 plane waves for each molecular orbital in the 20 a.u 
cubic cell. We performed our calculations without imposing periodic bound- 
ary conditions. Because we are dealing with charged isolated clusters with the 
strong dipoles and long-range Coulomb interactions, an aperiodic convolution 
method for solving Poisson's equations for free-space boundary conditions has 
been used [18]. We have also replaced hydrogen atoms by deuterium in our 
simulations in order to be able to use sufficiently large values for the time step 
{At — 7 a.u) and fictitious mass (// = 1100 a.u.) in Car - Parrinello dynamics. 



3 Results and discussion 

The acidity of a hydrated metal ion depends on the strength of the coordi- 
nation bond between cation and the oxygen atoms of the first hydration shell 
water molecules. If this coordination bond is strong, the OH bond in the sol- 
vating water molecules of the first coordination shell are strongly polarized. 
These polarized bond can dissociate in a polar solvent to form hydrolyzed 
species such as Al(H20)5(OH)^"'". This process is important to many aspects 
of aluminum solvation and is complicated by the many interactions that deter- 
mine the structure of the coordination shell. In the following, we will show that 
the polarization of the solvating molecules is a function of the coordination 
number of the first solvation shell. This leads to a cooperative relation be- 
tween mobility of the protons in the first solvation shell and the coordination 
structure. 
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3.1 Structure of aqueous Al^+ 



The formation of a stable coordination cell around Al^+ occurs on the subpi- 
cosecond time scale and hence well within the range of an ah initio molecular 
dynamics (AIMD) run. Our AIMD runs were initiated by placing one alu- 
minum ion in the middle of the cluster of 16 randomly oriented water molecules 
in a cubic simulation box with a length of 20 a.u. On the time scale of about 
0.25 ps, dramatic rearrangement of water cluster took place with the forma- 
tion of an almost perfect octahedra complex of A1(H20)6"'" surrounded by the 
10 remaining water molecules. This octahedral structure appears to be very 
stable and remained intact during the 1 ps run at ambient conditions. During 
this run no exchange of water molecules between first and second solvation 
shell was observed. 

To check the structural parameters of the hexa-hydrate complex as calculated 
by our method, we performed AIMD simulations on isolated A1(H20)q^ (by 
removing 10 remaining waters from the structure obtained in our first run) 
and also small clusters of [A1(H20)„]^"'" with n = 1,2,3,4,5. Since there are 
several possible low energy isomers of these clusters we have used simulating 
annealing, i.e., heating and cooling the system by reducing the kinetic energy, 
together with the steepest descent geometry optimizations to ensure that the 
lowest energy structures were obtained. In the present calculation we started 
from an initial temperature of 1000° K in the simulated annealing procedure. 
By performing short (on the time scale r = 10 fs) molecular dynamics runs 
followed by periodic quenching of the system, a rather fast effective cooling 
schedule was enforced. The geometrical parameters obtained in our calcula- 
tions (see Table 2) are in good agreement with the second order M0ller-Plesset 
perturbation theory (MP2) results [9]. 

For n = 1, 2, 3 planar structures were obtained. The n = 4 cluster was a tetra- 
hedral complex. At n = 5, the structure of A1(H20)5''' was a square-based 
pyramid (aluminum sits slightly above the center of the square formed by four 
waters with fifth water molecule along the line perpendicular to the square). 
This is different from Ref. [9], where A1(H20)5"'" was found to be a stable trigo- 
nal bipyramid (which can be viewed as distorted square-based pyramid) under 
constrained optimization. We tried to find a C2v bipyramid structure but 
the system always went to the square pyramid on annealing. The minimum en- 
ergy geometries of the trigonal bipyramid A1(H20)5^ and octahedral complex 
A1(H20)6"'" are given in the Fig. 1. 

Filling the first solvation shell of the aluminum cation by water molecules leads 
to gradual increase in the length of the coordination bond Al-0. At the same 
time the OH bond length of the waters in the first solvation shell approaches 
the value for free H2O at n = 6. The intramolecular angle < HOH for the 
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molecules in the first solvation shell is larger then in a free water molecule. 
As a function of coordination number the < HOH first decreases and then 
increases with the minimum at n = 3. This trend is also observed in MP2 
data but the numbers are smaller [9] . 

The cohesive energies per coordination bond for A1(H20)^"'" (n = 1 — 6) clus- 
ters have been calculated according to the formula 

AE^ = {E{A\'+{R20)^ - E{A\'+) - nE(H20))/n (1) 



using PBE96 functional (PBE96 results differed by a few percent from the 
magnitude of the LDA results). The cohesive energies per coordination bond 
A£^n as well as total cohesive energies AE^ x n for each of these clusters are 
reported in the Table 2 and are in close agreement with MP2 results [9] . 

To characterize the strength of polarization of the hydrated water molecules, in 
addition to the geometrical parameters of the optimized structures of [A1(H20)„]^"'" 
we also determined the approximate dipole moments of water molecules in 
these clusters. We found the centers (re) of the electronic density for the par- 
ticular water molecule according to the formulae 



where the integration is over a sphere of about 3 a.u. centered on the oxygen 
ions. We note that the electron density was well localized on the solvating 
waters. Using this procedure the dipole moment of isolated water molecule 
was found to be 1.838 D in good agreement with the experimental value of 
1.855 D [21]. The average dipoles of the single hydrated waters in [Al(H20)n]^''' 
clusters with n = 1, 2, 3, 4, 5, 6 decrease by almost 20% from 5.9 D to 5.02 D 
as a function of coordination. These numbers represent a substantial increase 
from the dipole moments of water molecule in vacuum and in the bulk liquid 
which are 1.855 D [21] and 2.6 D [22], respectively. 

The increase in the bond polarization as a function of hydration numbers plays 
an important role in the dynamical behavior of the systems as will be discussed 
below. To illustrate the reduction of energy in the transfer of a proton from 
water in the first solvation shell to vicinity of an oxygen in the second solvation 
shell in the presence of aqueous AP+ we carried out the following calculation. 
In the low energy structure of A1(H20)7'^, a proton was moved from the water 
in the first solvation shell to the nearest water molecule in the second solvation 
shell along hypothetical reaction path leading to the acid reaction 

Al3+(H20)6+H20(2nd solvation shell) — > A1^+(H20)50H-+H30+. (3) 
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The procedure was done at the constant 0-0 separation of 2.403 A. We kept 
0-H-O atoms colhnear while moving the proton between oxygen atoms in 
steps of 0.02 A. The potential energy calculated using the PBE96 functional 
along the path 0-H-O is plotted in Fig. 2 as a function of the asymmetric 
stretch coordinate 

^ _ koiul - \r 02h| / .V 

^ ~ ^/2 



where roiH and roiH are the distances between the transferring proton and 
oxygens of water molecules. The potential energy for the proton transfer pro- 
cess from the aqueous AP"^ to the near vicinity of the water molecule in the 
second solvation shell is represented by the broad potential well of 9 kcal/mol. 
We note that there is only a single minimum in the potential given in Fig. 2. To 
form a minimum corresponding to HaO^"*" species may require a more complex 
reaction coordinate [23]. This will be the subject of further research. A similar 
calculation for proton transfer in a water dimer was carried out. In this case 
the transfer of the proton from one water molecule to the other requires about 
50 kcal/mol. This dramatic reduction in the potential energy landscape due 
to the presence of aqueous Al^"*" contributes to high mobility of the protons in 
our molecular dynamics simulations of aluminum ion with n = 14 — 16 water 
molecules discussed below (see Section 3.3). 



3.2 Structure of aluminate ion, Al(OH) 



Since the LDA geometry of aqueous aluminum ion, A1(H20)q^ compared well 
with MP2 data and experimental X-ray diffraction results, we initiated a study 
of the formation of the most important hydrolysis product, A1(0H)4 . To do 
this we started an AIMD run with the equilibrated A1(H20)6''' structure but 
with 4 protons removed. The resulting deprotonatcd hexa-hydratc complex, 
A1(0H)4 (H20)2, after 0.1 ps readily evolved towards the equilibrium tetra- 
hedral structure of aluminate ion, Al(OH)j, see Fig. 3. The two remaining 
water molecules were forced out of the first solvation shell and formed hydro- 
gen bonds to the hydroxyl groups coordinated to the aluminate anion. 

Since these simulations were done at the finite ionic temperature, it allowed us 
to test the thermal stability of this process. The resulting negatively charged 
fourfold coordinated structure has a negative highest occupied eigenvalue of 
Kohn-Sham orbital of £ = —2.07 eV. The tetrahedral structure of A1(0H)J 
remained stable with no exchange in the following 0.5 ps of the AIMD run. 
Extensive MP2/6-311G ** and B3LYP/6-311G ** calculations of Al(OH)j 
gave the average Al-0 distance as 1.792 A [24] which compares well with 
our LDA result, 1.750 A. The Al-0 bonds are strengthened in tetrahedral 
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A1(0H)J structure as compared to the Al(H20)e''' octahedral complex, where 
the Al-0 distance is 9% larger. As was discussed in Section 1, the change in 
coordination from six in A1(H20)|^ to four in Al(OH)j is a unique feature of 
aluminum hydrolysis. 



3.3 Beyond first solvation shell 

The interaction between first and second solvation shells is an important issue 
in theoretical [8,9] and experimental studies [25,26] of the solvation of highly 
charged ions in water. As a first step to the understanding of these effects in 
bulk water, we report here results for the solvation of aluminum ion in the 
clusters of n = 14 — 16 waters. 

The presence of hydrolysis (e.g. acid reaction of the form of Eq. 3) in the 
AP+ water system is well documented [2]. However, our low temperature sim- 
ulations with AP+ in the middle of cluster of water molecules did not show 
any proton mobility on the time scale of the simulations. Since the hydrol- 
ysis constants in the aqueous aluminum ion solutions are known to be four 
orders of magnitude higher at the temperature of 800° K then at room tem- 
perature [27], wc also performed simulations with the Al"^+ inside cluster of 
14 water molecules in this elevated temperature regime. Hydrolysis did occur 
in this heated system on the subpicosecond time scale. As the coordination 
bond of one of the six hydrated waters (with purple colored oxygen) length- 
ened two protons (green) on the remaining waters of the first solvation shell 
shuttled between hydrated waters and the incomplete second solvation shell. 
Snapshots showing the shuttling protons coordinated to second solvation shell 
water molecules are presented in Fig. 4. 

In our calculations of the polarization of the OH bonds in the first coordination 
shell we have shown that the decrease in repulsive interactions in the solvation 
shell for low coordination led to shorter Al-0 bond lengths and corresponding 
increase in polarization of the OH bonds. In a polarizable medium (e.g. water) 
this increase in polarization should lead to increased mobility of the protons. 
Evidence for this from our simulations is given in Fig. 5 where we show a 
cluster structure in which the Al"^+ ion is coordinated to only 3 waters. The 
extensive migration of the protons even at room temperature shown in this 
figure support the above conjecture. Neighboring waters of AP"*" lose their 
protons shortly (~ 50 fs) after the cation has been artificially moved into 
the 3-water coordination structure on the surface of the water cluster. After 
equilibration at room temperature, one proton was found to shuttle between 
hydrated water and the water in the second solvation shell, while two other 
protons left first solvation shell and became engaged in the hydrogen bonding 
network of the solvent. These three protons which showed high mobility are 
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colored green in Fig. 5. Although this "surface" structure of aqueous aluminum 
is an artificial construction, it suggests the crucial role of the coordination of 
AP"*" in hydrolysis. 
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Table 1 

Intramolecular structural parameters for water monomer and hydrogen bond pa- 
rameters for water dimer in comparison with known LDA results and experimental 
data (in A and deg). 



calculated results Ref. [17] experiment 

water monomer 

roH 0.953 0.978 0.957" 

<HOH 104.8 104.4 104.5" 
water dimer 

roo 2.713 2.710 2.98 ^ 

< O - II.. .0 170.8 171 171 ±10 '' 

") Prom Ref. [19]. Prom Ref. [20]. 
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Table 2 

Global minimum geometries (in A and deg) and average dipole moments of water 
molecules (in Debye) of [Al(H20)n]^''" complexes with n= 1,2,3,4,5,6, obtained 
with the LDA density functional as described in the text. Total and per coordi- 
nation bond cohesive energies (in kcal/mol), AE'n x n and A.E^ respectively, were 
calculated using PBE96 functional. For comparison, all electron MP2/cc-pwCVTZ 
results for n = 1,6 and MP2/cc-pVDZ results for n = 5 from Ref.[9] are reported 
in parentheses. For the n = 6 the experimental X-ray diffraction result of Ref.[4] is 
given. 





n = 1 


n = 2 


n = 3 


n = 4 


n = 5 


n = 6 


TAIO 


1.724 


1.751 


1.751 


1.815 


1.823 1.870 1.880 


1.920 




(1.723) 








(1.895) (1.901) (1.949) 


(1.911) 














l.(S9 ((\\'i")erimoiii) 


roH 


1.031 


1.011 


0.999 


0.989 


0.980 


0.976 




(1.018) 








(0.984) 


(0.972) 


< HOH 


106.9 


106.5 


106.2 


106.7 


107.1 


107.6 




(106.54) 








(106.01) 


(106.56) 


dipole 


5.9 


5.51 


5.26 


5.06 


5.02 


5.02 




-203.7 


-184.3 


-166.45 


-148.16 


-131.25 


-118.36 


X n 


-203.7 


-368.6 


-499.35 


-592.64 


-656.27 


-710.16 




(-201.3) 








(-664.2) 


(-723.7) 
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Fig. 1. The minimum energy geometries of a) A1(H20)5 and b) Al(H20)g 



Fig. 2. Potential energy (in kcal/mol) for the proton moved between two water 
molecules in steps of 0.02 A along a linear path 0-H-O in Al(H20)y''" cluster. 0-0 
separation is 2.403 A. 0-H-O atoms have been kept collinear. Q is the asymmetric 
stretch coordinate (see text). 



Fig. 3. Snapshots from ab initio molecular dynamics simulation starting from the 

dcprotonatcd Al(OH)^ (H20)2 cluster, (a) Octahedral initial configuration; (b) after 
0.025 ps; (c) after 0.05 ps; and (d) the final tetrahedral structure of A1(0H)4 at 
0.3 ps with two remaining water molecules forced out of the first solvation shell. 
Hydrogen bonds are illustrated by dashed lines. 



Fig. 4. Snapshot of the hydrated Al(H20)g^ at a temperature of about 800° K with 
the protons (green) shuttling between hydrated waters. 



Fig. 5. Snapshot of the AIMD simulation of the artificial "surface" structure with 
a low coordinated aluminum ion on the surface of the water cluster and extensive 
hydrolysis. Protons which showed high mobility are colored green. 
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Fig. 1 b) 
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